Modeling the nonlinear hysteresis response of reservoir media

ABSTRACT

A method of modeling a nonlinear hysteresis response of reservoir media involves obtaining measurements of core samples from the reservoir to generate a multitude of bounding curves representing primary drainage and primary imbibition characteristics of a reservoir, generating, using a processor, a plurality of constitutive equations based on the multitude of bounding curves and a multitude of material functions, the multitude of material functions determined based on material properties associated with the reservoir, generating, using the processor, an output in a form of capillary-pressure scanning curves based on the multitude of constitutive equations, the capillary-pressure scanning curves representing drainage and imbibition behaviors of the reservoir during an oilfield operation, and storing the output.

CROSS-REFERENCE TO RELATED APPLICATIONS

This Application claims priority under 35 U.S.C. §119(e) to ProvisionalPatent Application No. 61/077,440, filed Jul. 1, 2008, entitled“GENERALIZED APPROACH FOR MODELING THE NONLINEAR HYSTERESIS RESPONSE OFRESERVOIR MEDIA,” which is incorporated herein by reference in itsentirety.

BACKGROUND

One of the fundamental reasons for the hysteretic nonlinear behavior ofporous reservoir media is that heterogeneous or damaged materialscontain an enormous number of mesoscopic features such as microcracksand macrocracks, joints, and grain to grain contacts containingmultiphase liquids. Each of these mesoscopic units exhibits a hystereticbehavior which dominates the macroscopic reservoir response. Fluids inthe reservoir may significantly influence the quasistatic and dynamicresponse in porous materials of the reservoir media due to theactivation of internal molecular forces. Hysteresis can affect oilrecovery efficiency because of the induced relative permeabilityreduction and the associated trapping of oil.

SUMMARY

In general, aspects of the invention are directed to a method ofmodeling a nonlinear hysteresis response of reservoir media involvesobtaining measurements of core samples from the reservoir to generate amultitude of bounding curves representing primary drainage and primaryimbibition characteristics of a reservoir generating using a processor aplurality of constitutive equations based on the multitude of boundingcurves and a multitude of material functions, the multitude of materialfunctions determined based on material properties associated with thereservoir, generating, using the processor, an output in a form ofcapillary-pressure scanning curves based on the multitude ofconstitutive equations, the capillary-pressure scanning curvesrepresenting drainage and imbibition behaviors of the reservoir duringan oilfield operation, and storing the output.

Other aspects of modeling the nonlinear hysteresis response of reservoirmedia will be apparent from the following description and the appendedclaims.

BRIEF DESCRIPTION OF THE DRAWINGS

So that the above recited features of modeling a nonlinear hysteresisresponse of reservoir media can be understood in detail, a moreparticular description, briefly summarized above, may be had byreference to the embodiments thereof that are illustrated in theappended drawings. It is to be noted, however, that the appendeddrawings illustrate typical embodiments of and are therefore not to beconsidered limiting of its scope, for modeling a nonlinear hysteresisresponse of reservoir media may admit to other equally effectiveembodiments.

FIGS. 1.1-1.4 depict a schematic view of an oilfield having subterraneanstructures containing reservoirs therein, in which embodiments ofmodeling a nonlinear hysteresis response of reservoir media can beimplemented.

FIGS. 2.1-2.4 depict graphical depictions of data collected by the toolsof FIGS. 1.1-1.4, respectively.

FIG. 3 depicts a schematic view, partially in cross section of anoilfield having a plurality of data acquisition tools positioned atvarious locations along the oilfield for collecting data from thesubterranean formations.

FIG. 4 shows an example schematic view of an oilfield having a pluralityof wellsites for producing hydrocarbons from the subterranean formation.

FIG. 5 shows an example schematic view of modeling a nonlinearhysteresis response of reservoir media in accordance with one or moreembodiments.

FIG. 6 is a flow chart of a method for modeling a nonlinear hysteresisresponse of reservoir media in accordance with one or more embodiments.

FIG. 7 depicts a computer system in accordance with one or moreembodiments of modeling a nonlinear hysteresis response of reservoirmedia.

DETAILED DESCRIPTION

Embodiments of modeling a nonlinear hysteresis response of reservoirmedia are shown in the above-identified figures and described in detailbelow. In describing the embodiments, like or identical referencenumerals are used to identify common or similar elements. The figuresare not necessarily to scale and certain features and certain views ofthe figures may be shown exaggerated in scale or in schematic in theinterest of clarity and conciseness.

In general, modeling a nonlinear hysteresis response of reservoir mediauses a set of constitutive equations based on measured bounding curvesto model hysteresis behaviors in material properties such as capillarypressure and relative permeability. In particular, the constitutiveequations may be configured based on material specific functions(referred to as material functions) to match laboratory measurements(e.g., of core samples from the reservoir) such that the modelingaccuracy can be fine-tuned.

FIGS. 1.1-1.4 show a schematic view of an oilfield having subterraneanstructures containing reservoirs therein, with various oilfieldoperations being performed on the oilfield.

FIG. 1.1 depicts a survey operation being performed to generate aseismic data output record (124) using recording truck computer (122.1)on a seismic recording truck (106.1) to receive, via geophone-receivers(118), data (120) of sound vibration(s) (112) that reflect off horizons(114) in an earth formation (116) from an acoustic source (110).

FIG. 1.2 depicts a drilling operation being performed by a drilling tool(106.2) suspended by a rig (128) and advanced into the subterraneanformation (102) to form a wellbore (136) for reaching the reservoir(104). Drilling mud is circulated through the drilling tool (106.2) viaa flow line (132) back to a mud pit (130) on the surface. The drillingtool may be adapted for measuring downhole properties such as adaptedfor taking a core sample (133). A surface unit (134) with a transceiver(137) collects data output (135) generated during the drilling operationand allows communications between various portions of the oilfield (100)or other locations.

FIG. 1.3 depicts a wireline operation and includes all the elementsdepicted in FIG. 1.2 except that the drilling tool (106.2) issubstituted with a wireline tool (106.3) adapted for performing welllogs, downhole tests, collecting samples, and/or performing a seismicsurvey operation based on an explosive or acoustic energy source (144)in which case the wireline tool (106.3) may provide data output (135) tothe surface unit (134).

FIG. 1.4 depicts a production operation being performed by a productiontool (106.4) deployed from a production unit or Christmas tree (129) andinto the completed wellbore (136) of FIG. 1.3 for drawing fluid from thedownhole reservoirs (104) into surface facilities (142) via a gatheringnetwork (146). Sensors (S) positioned about the oilfield (100) areoperatively connected to a surface unit (134) with a transceiver (137)for collecting data (135), for example, reservoir data, wellbore data,surface data and/or process data.

While one wellsite is shown, it will be appreciated that the oilfield(100) may cover a portion of land that hosts one or more wellsites.Part, or all, of the oilfield may be on land and/or sea. Also, theoilfield operations depicted in FIGS. 1.1-1.4 may be performed with anycombination of one or more oilfields, one or more processing facilitiesand one or more wellsites.

FIGS. 2.1-2.4 are graphical depictions of data collected by the tools ofFIGS. 1.1-1.4, respectively. FIG. 2.1 depicts a seismic trace (202) ofthe subterranean formation (102) of FIG. 1.1 taken by survey tool(106.1). FIG. 2.2 depicts a core sample (133) taken by the logging tool(106.2) of FIG. 1.2. FIG. 2.3 depicts a well log (204) of thesubterranean formation (102) taken by the wireline tool (106.3) of FIG.1.3. FIG. 2.4 depicts a production decline curve (206) of fluid flowingthrough the subterranean formation (102) taken by the production tool(106.4) of FIG. 1.4.

FIG. 3 is a schematic view, partially in cross section of an oilfield(300) having data acquisition tools (302.1), (302.2), (302.3), and(302.4) positioned at various locations along the oilfield forcollecting data of a subterranean formation (304). The data acquisitiontools (302.1-302.4) may be the same as data acquisition tools(106.1-106.4) of FIG. 1, respectively. As shown, the data acquisitiontools (302.1-302.4) generate data plots or measurements (308.1-308.4),respectively.

Data plots (308.1-308.3) are examples of static data plots that may begenerated by the data acquisition tools (302.1-302.4), respectively.Static data plot (308.1) is a seismic two-way response time and may bethe same as the seismic trace (202) of FIG. 2.1. Static plot (308.2) iscore sample data measured from a core sample of the formation (304),similar to the core sample (133) of FIG. 2.2. Static data plot (308.3)is a logging trace, similar to the well log (204) of FIG. 2.3. Data plot(308.4) is a dynamic data plot of the fluid flow rate over time, similarto the graph (206) of FIG. 2.4. Other data may also be collected, suchas historical data, user inputs, economic information, other measurementdata, and other parameters of interest.

The subterranean formation (304) has a plurality of geologicalstructures (306.1-306.4). As shown, the formation has a sandstone layer(306.1), a limestone layer (306.2), a shale layer (306.3), and a sandlayer (306.4). A fault line (307) extends through the formation. Thestatic data acquisition tools are adapted to measure the formation anddetect the characteristics of the geological structures of theformation.

While a specific subterranean formation (304) with specific geologicalstructures are depicted, it will be appreciated that the formation maycontain a variety of geological structures. Fluid may also be present invarious portions of the formation. Each of the measurement devices maybe used to measure properties of the formation and/or its underlyingstructures. While each acquisition tool is shown as being in specificlocations along the formation, it will be appreciated that one or moretypes of measurement may be taken at one or more location across one ormore oilfields or other locations for comparison and/or analysis.

The data collected from various sources, such as the data acquisitiontools of FIG. 3, may then be evaluated. Typically, seismic datadisplayed in the static data plot (308.1) from the data acquisition tool(302.1) is used by a geophysicist to determine characteristics of thesubterranean formation (304). Core data shown in static plot (308.2)and/or log data from the well log (308.3) is typically used by ageologist to determine various characteristics of the geologicalstructures of the subterranean formation (304). Production data from theproduction graph (308.4) is typically used by the reservoir engineer todetermine fluid flow reservoir characteristics.

FIG. 4 shows an oilfield (400) for performing production operations. Asshown, the oilfield has a plurality of wellsites (402) operativelyconnected to a central processing facility (454). The oilfieldconfiguration of FIG. 4 is not intended to limit the scope of theinvention. Part or all of the oilfield may be on land and/or sea. Also,while a single oilfield with a single processing facility and aplurality of wellsites is depicted, any combination of one or moreoilfields, one or more processing facilities and one or more wellsitesmay be present.

Each wellsite (402) has equipment that forms a wellbore (436) into theearth. The wellbores extend through subterranean formations (406)including reservoirs (404). These reservoirs (404) contain fluids, suchas hydrocarbons. The wellsites draw fluid from the reservoirs and passthem to the processing facilities via surface networks (444). Thesurface networks (444) have tubing and control mechanisms forcontrolling the flow of fluids from the wellsite to the processingfacility (454).

Generally speaking, for an oil reservoir, capillary pressure typicallyexhibits four distinct characteristics upon saturation. For example, interms of water saturation for a water-wet porous medium where drainageand imbibition denote decreasing and increasing water saturation,respectively, these characteristics are:

1) Primary drainage: occurs over a geological timescale when oil entersa porous medium previously saturated completely with water (e.g., watersaturation equals 1), resulting in displacement of water and reducedwater saturation resulting in the formation of an oil reservoir.

2) Secondary imbibition: typically occurs when water is injected intothe reservoir to displace oil, resulting in increased water saturation.

3) Secondary drainage: typically occurs when gas is injected into thereservoir, often as part of a water-alternating-gas (WAG) improved oilrecovery strategy. The injection operations (e.g., water injectionand/or WAG improved oil recovery) may be the same as those described inreference to FIG. 1.4 above.

Primary imbibition: occurs when water displaces oil from a porous mediumpreviously saturated completely with oil (e.g., water saturation equals0). Primary imbibition and other characteristics of capillary pressuredescribed above can be measured in the laboratory by performingexperiments on core samples taken from the reservoir, such as the coresamples described in reference to FIGS. 1.2 and 2.2 above.

For porous reservoir media, the nature of hysteresis remainsphenomenological and implies the existence of hysteretic mesoscopicunits (HMU). In this simplified phenomenological hysteresis model, aporous reservoir is represented by an assemblage of capillary elementsthat behave hysteretically as a function of the fluid saturation. Asingle HMU may exist in one of three states: imbibition, drainage, andintermediate scan states. The behavior of an HMU is such that the HMU isinitially located on the drainage curve with drainage capillary pressurep_(c) ^(d), but moves to attain the imbibition capillary pressure p_(c)^(i) by traversing a scanning curve as the saturation increases. Uponreaching the imbibition curve it remains on the curve as the saturationcontinues to increase. When the saturation starts to decrease, theelement once again traverses a scanning curve, which can be differentfrom an imbibition scanning curve, returning to its original drainagecurve. A large number of HMU elements with differing (p_(c) ^(d), p_(c)^(i)) parameters provides a signature for a heterogeneous porousreservoir, which depends on the state of saturation and on the degree ofdamage to the material.

As described above, modeling a nonlinear hysteresis response ofreservoir media may be used in Reservoir Simulation Modeling. In one ormore embodiments, a region (e.g., a HMU) of the reservoir may berepresented by one or more simulation grid blocks in the reservoirsimulation model, where capillary pressure p_(c) of each grid block maybe calculated based on p_(c) ^(d) and p_(c) ^(i) of the regioncontaining the grid block. Additional details related to modeling thecharacteristics of the capillary pressure are provided below.

FIG. 5 shows an example schematic view of modeling a nonlinearhysteresis response of reservoir media, which may be performed in anoilfield, such as the oilfield (300) depicted in FIG. 3 above.Specifically, FIG. 5 shows a schematic representation of the generalizedalgorithm for calculating capillary pressure characteristics as afunction of water saturation in a reservoir of the oilfield (300). Inone or more embodiments, the calculation may be based on a generalizedphenomenological model to describe the hysteretic nonlinear response ofcapillary pressure and relative permeability. For example, the activehysteresis may be modeled by the solution of either a differential or anintegral equation.

As shown in FIG. 5, computer program (504) receives input data (502) togenerate output (506). In one or more embodiments, the computer program(504) is configured to model characteristics of a reservoir, forexample, using an analytical model or a numerical method to solve thedifferential or integral equation describing the active hysteresis.Additional details related to the analytical model or numerical methodused in modeling a nonlinear hysteresis response of reservoir media areprovided below.

As shown in FIG. 5, the characteristics of the reservoir may include arelationship between the capillary pressure and fluid saturation. In oneor more embodiments, the computer program (504) receives an initialvalue of the capillary pressure/fluid saturation (i.e., an initial state(502.5) α₀) and generates the output (506) representing subsequentbehavior of the capillary pressure as a function of fluid saturation(e.g., water saturation S_(w)) when the reservoir condition shifts fromthe initial state (502.5) α₀; for example, in an injection operation. Inthe example shown in FIG. 5, the capillary pressure is denoted as p_(c)along the vertical axis while the fluid saturation, specifically watersaturation in this example, is denoted as S_(w) along the horizontalaxis. In addition, the primary drainage characteristics of the capillarypressure as a function of fluid saturation is shown as a primarydrainage curve (504.1) starting from S_(w)=1 following decreasing watersaturation to end at S_(w)=S_(wi) where capillary pressure reaches amaximum P_(c) ^(max), the primary imbibition characteristics of thecapillary pressure as a function of fluid saturation is shown as aprimary imbibition curve (504.2) starting from S_(w)=0 followingincreasing water saturation to end at S_(w)=S_(wi)−S_(or) wherecapillary pressure reaches a minimum p_(c) ^(min). The secondarydrainage characteristics of the capillary pressure as a function offluid saturation is shown as a secondary drainage curve (504.3) startingfrom S_(w)=S_(wi)−S_(or) with capillary pressure at the minimum p_(c)^(min) following decreasing water saturation to end at S_(w)=S_(wi) withcapillary pressure at the maximum p_(c) ^(max). The secondary imbibitioncharacteristics of the capillary pressure as a function of fluidsaturation is shown as a secondary imbibition curve (504.4) startingfrom S_(w)=S_(wi) with capillary pressure at the maximum p_(c) ^(max)following increasing water saturation to end at S_(w)=S_(wi)−S_(or) withcapillary pressure at the minimum p_(c) ^(min). The arrows on thesecondary drainage curve (504.3) and the secondary imbibition curve(504.4) indicate the direction of decreasing fluid saturation andincreasing fluid saturation respectively.

In one or more embodiments, portions or all of the curves(504.1)-(504.4) may be preconfigured in the computer program (504). Inone example, the preconfigured curves are based on a pre-definedanalytical model and may not be user adjustable. In another example, thepreconfigured curves may be user adjustable based on at least a portionof the input data (502.1)-(502.4). More details related to determiningthe preconfigured curves based on user input data are provided below.

Further as shown in FIG. 5, an initial state (e.g., initial values ofthe capillary pressure and fluid saturation) of a grid block in thereservoir model may be represented by point A, which corresponds to aninitial value of a hysteretic function parameter α₀ (502.5). Additionaldetails related to the hysteretic function parameter are provided below.In a water or WAG injection operation, the state of the grid block maychange in time from point A as the capillary pressure and fluidsaturation change due to the injection operations. For example, thestate of the grid block may traverse a curve from point A toward point Eon the secondary imbibition curve (504.4). Based on the specificinjection schedule, the state of the grid block may reach a reversalpoint B and change course to traverse another curve from point B towardpoint F on the primary drainage curve (504.1). These curves traversed bythe state of the grid block as the reservoir condition changes arereferred to as scanning curves. At any point of time, the state of thegrid block may be in an intermediate scan state (e.g., points A, B, andC) when it is not on any of the primary or secondary drainage/imbibitioncurves (504.1)-(504.4).

In one or more embodiments, constitutive equations may be used to modelthe characteristics of the capillary pressure described above. Aconstitutive equation is a relation between two physical quantities thatis specific to a material or substance and approximates the response ofthat material to external forces. Some constitutive equations arephenomenological while others are derived from first principles (i.e.,basic assumptions). A common constitutive equation is expressed as asimple proportionality using a parameter taken to be a property of thematerial. One such example of a constitutive equation is Hooke's law.The constitutive equation may be combined with other equations governingphysical laws (e.g., Darcy's law, various conservation laws, etc.) tosolve physical problems.

In one or more embodiments, the constitutive equations for capillarypressure when bounded between the primary drainage (p^(d) _(c)) andimbibition (p^(i) _(c)) curves can be written as:

$\begin{matrix}{{p_{c} = {{\left( {1 - \alpha} \right)p_{c}^{d}} + {\alpha \; p_{c}^{i}}}},{\alpha \in \left\lbrack {0,1} \right\rbrack}} & (1) \\{{\frac{\alpha}{S} = {{\eta (\alpha)}\left\lbrack {{{{sgn}\left( \overset{.}{S} \right)}\left( {{\chi\left( {S,\overset{.}{S}} \right)} - \alpha} \right)} + {\gamma\left( {S,\overset{.}{S}} \right)}} \right\rbrack}},{\alpha \in \left\lbrack {0,1} \right\rbrack}} & (2)\end{matrix}$

where α is the hysteretic function parameter, S is the fluid saturation,sgn( ) is the sign function, and {dot over (S)}=dS/dt is the rate of thefluid saturation. Further, α, p_(c), p^(d) _(c), and p^(i) _(c) arefunctions of S but are abbreviated for clarity while χ (S,{dot over(S)}), γ(S,{dot over (S)}) and η(α) are material functions that may bespecified based on material properties to model phenomenologicalbehaviors of the reservoir. The primary drainage (p^(d) _(c)) andimbibition (p^(i) _(c)) curves are referred to as bounding curves forthe region represented by equation (1).

In one or more embodiments, the computer program (504) is configured touse constitutive equations (e.g., (1) and (2)) to model the secondarydrainage/imbibition curves and the scanning curves of the reservoir. Inone or more embodiments, the material functions χ(S,{dot over (S)}),γ(S,{dot over (S)}), η(α) and the parameters used in the model may befine-tuned to match different experimental data or can be used ashistory matching parameters.

In one or more embodiments, the material functions may be specified asbelow:

$\begin{matrix}{{\eta (\alpha)} = \left\{ \begin{matrix}{{\alpha \; \psi_{0}},} & {{{{if}\mspace{14mu} \frac{S}{t}} < {0 - {drainage}}},} \\{{\left( {1 - \alpha} \right)\psi_{0}},} & {{{if}\mspace{14mu} \frac{S}{t}} > {0 - {imbibition}}}\end{matrix} \right.} & (3) \\{{\chi\left( {S,\overset{.}{S}} \right)} = \left\{ {\begin{matrix}{{\chi_{1}(S)},} & {{{if}\mspace{14mu} \frac{S}{t}} < 0} \\{{\chi_{2}(S)},} & {{{if}\mspace{14mu} \frac{S}{t}} > 0}\end{matrix},{{\gamma\left( {S,\overset{.}{S}} \right)} = \left\{ \begin{matrix}{{\gamma_{1}(S)},} & {{{if}\mspace{14mu} \frac{S}{t}} < 0} \\{{\gamma_{2}(S)},} & {{{if}\mspace{14mu} \frac{S}{t}} > 0}\end{matrix} \right.}} \right.} & (4)\end{matrix}$

where ψ₀ is a material specific constant (e.g., 0.25).

Equations (1-4) describe the subsequent transition between imbibitionand drainage curves. Material functions may be specified such that thehysteresis models obey the symmetry restriction with respect to thetransformation (1−α)→α. For these types of hysteresis models, thefollowing equality restriction is obtained:

χ₂(S)+χ₁(S)≡1+γ₁(S)−γ₂(S)  (5)

In one or more embodiments, the material function may be furtherspecified as γ₁(S)=Kγ₂ (S); therefore, equation (4) can be rewritten as:

χ₂(S)+χ₁(S)≡1+(K−1)γ₂(S), γ₁(S)=Kγ ₂(S)  (6)

Finally, the equations (1-4) and (6) form the set of constitutiveequations for modeling a nonlinear hysteresis response of reservoirmedia. The set of constitutive equations are based on predefinedbounding curves, for example the primary drainage (p^(d) _(c)) andimbibition (p^(i) _(c)) curves measured from core samples. In one ormore embodiments, the computer program (504) is configured to generatethe output (506) using these constitutive equations (1-4), (6) and theassociated predefined bounding curves.

In one or more embodiments, the differential equation (2) for the casewhen {dot over (S)}≠0 can be solved analytically and the model expressedin terms of the bounding and scanning curves.

In one embodiment, the material functions

${\chi\left( {S,\overset{.}{S}} \right)} = {H\left( {- \frac{S}{t}} \right)}$

(where H(•) is the Heaviside step function) and γ(S,{dot over (S)})≡0are considered:

$\begin{matrix}{\frac{\alpha}{S} = \left\{ {\begin{matrix}{{\left( {1 - \alpha} \right)^{2}\psi_{0}},} & {\frac{S}{t} > 0} \\{{\alpha^{2}\psi_{0}},} & {\frac{S}{t} < 0}\end{matrix},{S \in \left\lbrack {S_{f}^{\min},S_{f}^{\max}} \right\rbrack}} \right.} & (7)\end{matrix}$

where S_(f) ^(min) is the minimum fluid saturation and S_(f) ^(max) isthe maximum fluid saturation.

The differential equation (7) with the defined initial data can besolved analytically. For the case when {dot over (S)}>0, the analyticalsolution has the form:

$\begin{matrix}{{{\alpha (S)} = \frac{\left( {\frac{1}{S - S_{f}^{1} - E} - \frac{1}{E}} \right)}{\left( {\frac{1}{S_{f}^{2} - S_{f}^{1} - E} - \frac{1}{E}} \right)}},{{{where}\mspace{14mu} {\alpha \left( S_{1} \right)}} = 0},{{\alpha \left( S_{2} \right)} = 1}} & (8)\end{matrix}$

and

$E = \frac{1}{\psi_{0}}$

is the curvature parameter.

In another embodiment, the material functions are specified by χ(S,{dotover (S)})≡½ and γ(S,{dot over (S)})≡0 providing a trivial solution toequation (5) and equation (2) takes the form:

$\begin{matrix}{\frac{\alpha}{S} = \left\{ {\begin{matrix}{{\left( {1 - \alpha} \right)\left( {\frac{1}{2} - \alpha} \right)\psi_{0}},} & {\frac{S}{t} > 0} \\{{{- {\alpha \left( {\frac{1}{2} - \alpha} \right)}}\psi_{0}},} & {\frac{S}{t} < 0}\end{matrix},{S_{w} \in \begin{bmatrix}{S_{wi},} & S_{wmax}\end{bmatrix}}} \right.} & (9)\end{matrix}$

It should be noted that equation (9) is valid for

${\alpha \in \left\lbrack {0,\frac{1}{2}} \right\rbrack},$

which does not satisfy the condition that αε[0,1].

In yet another embodiment, the material functions χ(S,{dot over (S)}),γ(S,{dot over (S)}) are defined as:

$\begin{matrix}{{\chi\left( {S,\overset{.}{S}} \right)} = \left\{ {\begin{matrix}{\frac{K}{2},} & {{{if}\mspace{14mu} \frac{S}{t}} < 0} \\{\frac{K}{2},} & {{{if}\mspace{14mu} \frac{S}{t}} > 0}\end{matrix},{{\gamma\left( {S,\overset{.}{S}} \right)} = \left\{ \begin{matrix}{K_{1},} & {{{if}\mspace{14mu} \frac{S}{t}} < 0} \\{K_{2},} & {{{if}\mspace{14mu} \frac{S}{t}} > 0}\end{matrix} \right.}} \right.} & (10)\end{matrix}$

where K₁=K and K₂=1 are assumed. The final exponential hysteresis modeltakes the form:

$\begin{matrix}{\frac{\alpha}{S} = \left\{ {\begin{matrix}{{\left( {1 - \alpha} \right)\left( {\frac{K + 2}{2} - \alpha} \right)\psi_{0}},} & {\frac{S}{t} > 0} \\{{{\alpha \left( {\alpha + \frac{K}{2}} \right)}\psi_{0}},} & {\frac{S}{t} < 0}\end{matrix},{S_{w} \in \left\lbrack {S_{wi},S_{wmax}} \right\rbrack}} \right.} & (11)\end{matrix}$

where K>0 is the exponential hysteresis parameter and αε[0,1]. Theresulting differential equation (11) with the defined initial data canbe solved analytically, for example, for {dot over (S)}>0:

$\begin{matrix}{{{\alpha (S)} = \frac{{\frac{K + 2}{2}{\exp\left( {\frac{K}{2}\left\lbrack {\frac{S}{E} - \frac{S_{f}^{1}}{E} - {\frac{2}{K}{\ln\left( \frac{K + 2}{2} \right)}}} \right\rbrack} \right)}} - 1}{{\exp\left( {\frac{K}{2}\left\lbrack {\frac{S}{E} - \frac{S_{f}^{1}}{E} - {\frac{2}{K}{\ln\left( \frac{K + 2}{2} \right)}}} \right\rbrack} \right)} - 1}},} & (12)\end{matrix}$

where α(S_(f) ¹)=0 and

$E = \frac{1}{\psi_{0}}$

is the effective curvature. Equations (1) and (12) constitute theexponential hysteresis model.

In still another embodiment, or the trigonometric hysteresis model, thematerial functions χ(S,{dot over (S)}) and γ(S,{dot over (S)}) areconsidered as follows:

$\begin{matrix}{{\chi\left( {S,\overset{.}{S}} \right)} = \left\{ {{\begin{matrix}{{K\; {\sin^{2}\left( {2\pi \; S} \right)}},} & {{{if}\mspace{14mu} \frac{S}{t}} < 0} \\{{K\; {\cos^{2}\left( {2\pi \; S} \right)}},} & {{{{if}\mspace{14mu} \frac{S}{t}} > 0},}\end{matrix}{\gamma\left( {S,\overset{.}{S}} \right)}} = \left\{ \begin{matrix}{K_{1},} & {{{if}\mspace{14mu} \frac{S}{t}} < 0} \\{K_{2},} & {{{if}\mspace{14mu} \frac{S}{t}} > 0}\end{matrix} \right.} \right.} & (13)\end{matrix}$

where K₁=K and K₂=1 are assumed. In this case, equation (5) is satisfieddue to the fundamental trigonometric equality:

χ₂(S)+χ₁(S)=K cos²(2πS)+K sin²(2πS)≡K  (14)

Therefore, equation (2) takes the form:

$\begin{matrix}{\frac{\alpha}{S} = \left\{ {{\begin{matrix}{{\left( {1 - \alpha} \right)\left( {{\cos^{2}\left( {2\pi \; S} \right)} + 1 - \alpha} \right)\psi_{0}},} & {\frac{S}{t} > 0} \\{{{\alpha \left( {K + \alpha - {\sin^{2}\left( {2\pi \; S} \right)}} \right)}\psi_{0}},} & {{\frac{S}{t} < 0},}\end{matrix}S_{w}} \in \left\lbrack {S_{wi},S_{wmax}} \right\rbrack} \right.} & (15)\end{matrix}$

where K>0 is the exponential hysteresis parameter and αε[0,1]. Theresulting differential equation (15) with the defined initial data canbe solved analytically, for example, for {dot over (S)}<0, resulting ina differential equation in Bernoulli's form:

$\begin{matrix}{{\frac{\alpha}{S} = {{\psi_{0}\alpha^{2}} + {{\alpha \left( {K - {\sin^{2}\left( {2\pi \; S} \right)}} \right)}\psi_{0}}}},{S_{w} \in \left\lbrack {S_{wi},S_{wmax}} \right\rbrack},} & (16)\end{matrix}$

which can be solved analytically by Bernoulli's substitution

$\frac{1}{\alpha} = {v\text{:}}$

$\begin{matrix}{{{\frac{1}{\alpha} = {{C(S)}{\exp \left( {- {\omega (S)}} \right)}}},{{\omega (S)} = {{{\psi_{0}\left( {K + \frac{1}{2}} \right)}S} - {\frac{1}{8\pi}\psi_{0}{\sin \left( {4\pi \; S} \right)}}}}}{{{C(S)} = {{{- \psi_{0}}{\int{{\exp \left( {- {\omega (S)}} \right)}{S}}}} + C_{1}}},}} & (17)\end{matrix}$

Although certain examples described above apply to the case where {dotover (S)}>0, one skilled in the art will appreciate that thedifferential equations for the case where {dot over (S)}<0 may be solvedin similar manner and vice versa.

The constitutive equations (1-4) and (6) provide a general mathematicalformulation for the hysteresis model. Based on specific materialproperties of the reservoir being modeled, the user may select one ofthe analytical models described above that produces output (506)sufficiently consistent with actual measured data.

The analytical solution described above is efficient for computingresources and performance; however, in practice, some form of usercontrol over the nature of the hysteresis may be useful. Thedifferential form of the equations (1-4) and (6) allows a smoothtransition, as observed in experimental data, for both relativepermeability and capillary pressure from drainage-to-imbibition andimbibition-to-drainage states and requires minimum storage of parametersduring the simulation and, thus, is suitable for incorporation into anumerical model. An appropriate choice of material functions may bebased on tabular data obtained from laboratory measurements, for exampleof core samples described in reference to FIGS. 1.2 and 2.2 above.

In one or more embodiments, the computer program (504) may be configuredto determined the bounding curves and material functions of theconstitutive equations under user control, for example, based on thefollowing experimental data as input data (502) from the user:(S_(f))_(i=1 . . . N), (p_(prim) ^(d))_(i=1 . . . N), (p_(prim)^(i))_(i=1 . . . N), (p_(sec) ^(d))_(i=1, . . . N), and (p_(sec)^(i))_(i=1 . . . N), where (S_(f))_(i=1 . . . N) is the experimentaldata for fluid saturations, (p_(prim) ^(d))_(i=1 . . . N) is theexperimental primary drainage curve, (p_(prim) ^(i))_(i=1 . . . N) isthe experimental primary imbibition curve, (p_(sec)^(d))_(i=1, . . . , N) is the experimental secondary drainage curve,(p_(sec) ^(i))_(i=1 . . . N) is the experimental secondary imbibitioncurve and N is the number of experimental points as depicted in FIG. 5.

In one or more embodiments, the computer program (504) may be configuredto determine the bounding curves based on (S_(f))_(i=1 . . . N)(p_(prim) ^(d))_(i=1 . . . N), and (p_(prim) ^(i))_(i=1 . . . N), forexample, using interpolation or intrapolation techniques known in theart.

In one or more embodiments, the computer program (504) may be configuredto determine the material functions using a numerical method to achievesufficient consistency between the modeled results versus the measuredresults, (S_(f))_(i=1 . . . N), (p_(sec) ^(d))_(i=1, . . . , N), and(p_(sec) ^(i))_(i=1 . . . N). Additional details related to thenumerical method are provided below.

Applying the second order central difference scheme, the followingapproximation of derivatives (dα/dS_(f)) can be obtained:

$\begin{matrix}{{{\left( \frac{\alpha}{S_{f}} \right)_{i}^{drn} = \frac{\alpha_{i + 1} - \alpha_{i - 1}}{\left( S_{f} \right)_{i + 1} - \left( S_{f} \right)_{i - 1}}},{i = {1\mspace{14mu} \ldots \mspace{14mu} N}}}{{\text{-}{Secondary}\mspace{14mu} {drainage}\mspace{14mu} {curve}};}} & (18) \\{{{\left( \frac{\alpha}{S_{f}} \right)_{i}^{imb} = \frac{\alpha_{i + 1} - \alpha_{i - 1}}{\left( S_{f} \right)_{i + 1} - \left( S_{f} \right)_{i - 1}}},{{i = {1\mspace{14mu} \ldots \mspace{14mu} N}};}}{{\text{-}{Secondary}\mspace{14mu} {imbibition}\mspace{14mu} {curve}};}} & (19)\end{matrix}$

Applying constitutive equations (2-4) for the secondary drainage andimbibition curves and equation (13) for the material function γ(S,{dotover (S)}), the following system of equations can be written:

$\begin{matrix}{{{\chi_{1}\left( S_{f}^{i} \right)} = {\alpha_{i} - {\frac{1}{\alpha_{i}\psi_{0}}\left( \frac{\alpha}{S_{f}} \right)_{i}^{drn}} + K_{1}}},{{\text{-}{Drainage}\mspace{14mu} {scanning}\mspace{14mu} {curve}};}} & (20) \\{{{\chi_{2}\left( S_{f}^{i} \right)} = {\alpha_{i} + {\frac{1}{\left( {1 - \alpha_{i}} \right)\psi_{0}}\left( \frac{\alpha}{S_{f}} \right)_{i}^{imb}} - K_{2}}},{{\text{-}{Imbibition}\mspace{14mu} {scanning}\mspace{14mu} {curve}};}} & (21)\end{matrix}$

where K₁ and K₂ are user defined constants. This algorithm can benumerically implemented and generates capillary-pressure scanning curveswith provisions for hysteresis loops from any reversal point.

Although specific examples are described in the embodiments above, oneskilled in the art will appreciate that the computer program (504) isnot constrained to the use of only drainage and imbibition curves, whichare assumed to form a closed loop hysteresis.

FIG. 6 is a flow chart of a method for modeling a nonlinear hysteresisresponse of reservoir media in accordance with one or more embodiments.In one or more embodiments, one or more of the elements shown in FIG. 6may be omitted, repeated, and/or performed in a different order.Accordingly, embodiments of modeling a nonlinear hysteresis response ofreservoir media should not be considered limited to the specificarrangements of elements shown in FIG. 6.

In one or more embodiments, the method depicted in FIG. 6 may bepracticed based on the schematic view described with respect to FIG. 5above. Initially, bounding curves representing primary drainage andprimary imbibition characteristics of the reservoir may be obtained(Element 601). As described above, the primary drainage and primaryimbibition characteristics may include hysteresis characteristics. Inone or more embodiments, the bounding curves may include the primarydrainage curve and the primary imbibition curve as described withrespect to FIG. 5 above, which may be measured based on cores samplesfrom the reservoir. For example, the measurements may be represented as(S_(f))_(i=1 . . . N), (p_(prim) ^(d))_(i=1 . . . N), and (p_(prim)^(i))_(i=1 . . . N).

In Element 602, empirical secondary drainage and imbibition curves maybe measured based on core samples from the reservoir. For example, themeasurements may be represented as (S_(f))_(i=1 . . . N), (p_(sec)^(d))_(i=1, . . . , N), and (p_(sec) ^(i))_(i=1 . . . N) describedabove.

In Element 603, material functions (e.g., equations (3)-(6)) may bedetermined based on material properties of the reservoir. In one or moreembodiments, the material functions may be pre-determined (e.g.,equations (7), (10), or (13)) for multiple specific reservoir types,where the user may select appropriate material functions for use in thereservoir simulation modeling based on the type of reservoir beingmodeled. In one or more embodiments, the material functions may bedetermined using a numerical method by matching numerically generatedsecondary drainage and imbibition curves (e.g., equations (18) and (19))to the empirical secondary drainage and imbibition curves (e.g., asmeasured in Element 602 above).

In Element 604, constitutive equations (e.g., equations (1), (2)) may begenerated based on the bounding curves and the material functions. Inone or more embodiments, the constitutive equations may be solvedanalytically (e.g., equations (8), (9), (12), or (17)). In one or moreembodiments, the constitutive equations may be solved using a numericalmethod (e.g., equations (20), (21)).

In Element 605, an output in the form of capillary-pressure scanningcurves may be generated based on the plurality of constitutiveequations. For example the capillary-pressure scanning curves may berepresented by equations (8), (9), (12), (17), (20), or (21). Asdescribed with respect to FIG. 5 above, the capillary-pressure scanningcurves represent drainage and/or imbibition behaviors of the reservoirduring the oilfield operation. For example, a scanning curve may startfrom an initial state of the reservoir condition at the onset of aninjection operation and represent the traversal of the reservoircondition as the injection operation continues according to theinjection schedule or strategy. Optionally for the analytical solution,the adequacy of the selected material functions may be confirmed bycomparing the secondary drainage and imbibition curves derived from theanalytical model to the empirical secondary drainage and imbibitioncurves (e.g., as measured in Element 602 above).

Embodiments of the invention may be implemented on virtually any type ofcomputer regardless of the platform being used. For example, as shown inFIG. 7, a computer system (700) includes one or more processor(s) (702),associated memory (704) (e.g., random access memory (RAM), cache memory,flash memory, etc.), a storage device (706) (e.g., a hard disk, anoptical drive such as a compact disk drive or digital video disk (DVD)drive, a flash memory stick, etc.), and numerous other elements andfunctionalities typical of today's computers (not shown). The computer(700) may also include input means, such as a keyboard (708), a mouse(710), or a microphone (not shown). Further, the computer (700) mayinclude output means, such as a monitor (712) (e.g., a liquid crystaldisplay (LCD), a plasma display, or cathode ray tube (CRT) monitor).Those skilled in the art will appreciate that many different types ofcomputer systems exist, and the aforementioned input and output meansmay take other forms. Generally speaking, the computer system (700)includes at least the minimal processing, input, and/or output meansnecessary to practice embodiments of the invention.

Further, those skilled in the art will appreciate that one or moreelements of the aforementioned computer system (700) may be located at aremote location and connected to the other elements over a network.Further, embodiments of the invention may be implemented on adistributed system having a plurality of nodes, where each portion ofthe invention may be located on a different node within the distributedsystem. In one embodiment of the invention, the node corresponds to acomputer system. Alternatively the node may correspond to a processorwith associated physical memory. The node may alternatively correspondto a processor with shared memory and/or resources. Further, softwareinstructions to perform embodiments of the invention may be stored on acomputer readable medium such as a compact disc (CD), a diskette, atape, a file, or any other computer readable storage device.

One or more of the embodiments of modeling a nonlinear hysteresisresponse of reservoir media may include one or more of the following.One or more of the embodiments may provide a more consistentmathematical background to existing empirical solutions of thehysteresis behavior observed on variables, such as the capillarypressure. The mathematical framework used in one or more of theembodiments may allow for a better analysis of the hysteresis behavior,and potentially lead to more original solutions based on an appropriatechoice of the materials functions in the model. Further, in one or moreof the embodiments, the mathematical framework models the hysteresis inrelative permeabilities. In one or more of the embodiments, a correctmodeling of capillary pressure and relative permeability hysteresis isutilized in the assessment of the quality of a reservoir simulator and,as a consequence, in the quality of reservoir production management. Inone or more of the embodiments, the proposed model allows a smoothtransition of both relative permeability and capillary pressure fromdrainage-to-imbibition or imbibition-to-drainage states and requiresminimal storage of parameters during the simulation. One skilled in theart will appreciate that the characteristics described above are notexhaustive and the embodiments should not be limited by those provided.

It will be understood from the foregoing description that variousmodifications and changes may be made in the embodiments of modeling anonlinear hysteresis response of reservoir media without departing fromits true spirit. For example, the simulators couplings and arrangementof the system may be selected to achieve the desired simulation. Thesimulations may be repeated according to the various configurations, andthe results compared and/or analyzed.

This description is intended for purposes of illustration and should notbe construed in a limiting sense. The scope of this invention should bedetermined by the language of the claims that follow. The term“comprising” within the claims is intended to mean “including at least”such that the recited listing of elements in a claim are an open group.“A,” “an” and other singular terms are intended to include the pluralforms thereof unless specifically excluded.

1. A method of modeling a nonlinear hysteresis response of reservoirmedia comprising: obtaining measurements of core samples from thereservoir to generate a plurality of bounding curves representingprimary drainage and primary imbibition characteristics of a reservoir;generating, using a processor, a plurality of constitutive equationsbased on the plurality of bounding curves and a plurality of materialfunctions, the plurality of material functions determined based onmaterial properties associated with the reservoir; generating, using theprocessor, an output in a form of capillary-pressure scanning curvesbased on the plurality of constitutive equations, the capillary-pressurescanning curves representing drainage and imbibition behaviors of thereservoir during an oilfield operation; and storing the output.
 2. Themethod of claim 1, wherein the capillary-pressure scanning curves aregenerated using an analytical solution of the plurality of constitutiveequations.
 3. The method of claim 2, further comprising: identifying thematerial functions based on characteristics of the reservoir.
 4. Themethod of claim 1, wherein the capillary-pressure scanning curves aregenerated using a numerical solution of the plurality of constitutiveequations, wherein the material functions are determined by matchingnumerically generated secondary drainage and imbibition curves toempirical secondary drainage and imbibition curves measured from coresamples from the reservoir.
 5. The method of claim 1, wherein theplurality of constitutive equations comprises:p _(c)=(1−α)p _(c) ^(d) +αp _(c) ^(i),αγ[0,1],dα/dS=?η(α)[sgn({dot over (S)})(χ(S,{dot over (S)})−α)+γ(S,{dot over(S)})], αε[0,1] where p_(c) is capillary pressure, p_(c) ^(d) is aprimary drainage curve, p_(c) ^(i) is a primary imbibition curve, α is ahysteretic function parameter, S is fluid saturation, sgn( ) is signfunction, {dot over (S)}=dS/dt is a rate of the fluid saturation, andχ(S,{dot over (S)}), γ(S,{dot over (S)}) and η(α) are at least a portionof the plurality of material functions.
 6. The method of claim 5,wherein the plurality of material functions comprises:${\eta (\alpha)} = \left\{ {{\begin{matrix}{{\alpha\psi}_{0},} & {{{if}\mspace{14mu} \frac{S}{t}} < 0} & {{\text{-}{drainage}},} \\{{\left( {1 - \alpha} \right)\psi_{0}},} & {{{if}\mspace{14mu} \frac{S}{t}} > 0} & {{\text{-}{imbibition}},}\end{matrix}{\chi\left( {S,\overset{.}{S}} \right)}} = \left\{ {\begin{matrix}{{\chi_{1}(S)},} & {{{if}\mspace{14mu} \frac{S}{t}} < 0} \\{{\chi_{2}(S)},} & {{{{if}\mspace{14mu} \frac{S}{t}} > 0},}\end{matrix},{{\gamma\left( {S,\overset{.}{S}} \right)} = \left\{ \begin{matrix}{{\gamma_{1}(S)},} & {{{if}\mspace{14mu} \frac{S}{t}} < 0} \\{{\gamma_{2}(S)},} & {{{if}\mspace{14mu} \frac{S}{t}} > 0}\end{matrix} \right.}} \right.} \right.$ where ψ₀ is a material specificconstant and t represents time.
 7. The method of claim 6, wherein theplurality of material functions further comprises:χ₂(S)+χ₁(S)≡1+(K−1)γ₂(S), γ₁(S)=Kγ₂(S), where K is a constant.
 8. Acomputer readable medium for modeling a nonlinear hysteresis response ofreservoir media, the computer readable medium storing instructions whenexecuted by the processor comprising functionality for: obtaining aplurality of bounding curves representing primary drainage and primaryimbibition characteristics of a reservoir; generating, using aprocessor, a plurality of constitutive equations based on the pluralityof bounding curves and a plurality of material functions, the pluralityof material functions determined based on material properties associatedwith the reservoir; obtaining, using the processor, an analyticalsolution of the plurality of constitutive equations ascapillary-pressure scanning curves, the capillary-pressure scanningcurves representing drainage and imbibition behaviors of the reservoirduring an oilfield operation; and storing the output.
 9. The computerreadable medium of claim 8, wherein the plurality of bounding curves aregenerated based on measurements of core samples from the reservoir. 10.The computer readable medium of claim 8, wherein the material functionsare pre-determined based on characteristics of the reservoir.
 11. Thecomputer readable medium of claim 8, wherein the capillary-pressurescanning curves are generated using a numerical solution of theplurality of constitutive equations, wherein the material functions aredetermined by matching numerically generated secondary drainage andimbibition curves to empirical secondary drainage and imbibition curvesmeasured from core samples from the reservoir.
 12. The computer readablemedium of claim 8, wherein the plurality of constitutive equationscomprises:${p_{c} = {{\left( {1 - \alpha} \right)p_{c}^{d}} + {\alpha \; p_{c}^{i}}}},{\alpha \in \left\lbrack {0,1} \right\rbrack},{\frac{\alpha}{S} = {{\eta (\alpha)}\left\lbrack {{{{sgn}\left( \overset{.}{S} \right)}\left( {{\chi\left( {S,\overset{.}{S}} \right)} - \alpha} \right)} + {\gamma\left( {S,\overset{.}{S}} \right)}} \right\rbrack}},{\alpha \in \left\lbrack {0,1} \right\rbrack}$where p_(c) is capillary pressure, p_(c) ^(d) is a primary drainagecurve, p_(c) ^(i) is a primary imbibition curve, α is a hystereticfunction parameter, S is fluid saturation, sgn( ) is sign function, {dotover (S)}=dS/dt is a rate of the fluid saturation, and χ(S,{dot over(S)}), γ(S,{dot over (S)}) and η(α) are at least a portion of theplurality of material functions.
 13. The computer readable medium ofclaim 12, wherein the plurality of material functions comprises:${\eta (\alpha)} = \left\{ {{\begin{matrix}{{\alpha\psi}_{0},} & {{{if}\mspace{14mu} \frac{S}{t}} < 0} & {{\text{-}{drainage}},} \\{{\left( {1 - \alpha} \right)\psi_{0}},} & {{{if}\mspace{14mu} \frac{S}{t}} > 0} & {{\text{-}{imbibition}},}\end{matrix}{\chi\left( {S,\overset{.}{S}} \right)}} = \left\{ {\begin{matrix}{{\chi_{1}(S)},} & {{{if}\mspace{14mu} \frac{S}{t}} < 0} \\{{\chi_{2}(S)},} & {{{{if}\mspace{14mu} \frac{S}{t}} > 0},}\end{matrix},{{\gamma\left( {S,\overset{.}{S}} \right)} = \left\{ \begin{matrix}{{\gamma_{1}(S)},} & {{{if}\mspace{14mu} \frac{S}{t}} < 0} \\{{\gamma_{2}(S)},} & {{{if}\mspace{14mu} \frac{S}{t}} > 0}\end{matrix} \right.}} \right.} \right.$ where ψ₀ is a material specificconstant and t represents time.
 14. The computer readable medium ofclaim 13, wherein the plurality of material functions further comprises:χ₂(S)+χ₁(S)≡1+(K−1)γ₂(S), γ₂(S)=Kγ ₂(S), where K is a constant.
 15. Acomputer system for modeling a nonlinear hysteresis response ofreservoir media, the system comprising: a processor, memory storinginstructions when executed by the processor comprising functionalityfor: obtaining a plurality of bounding curves representing primarydrainage and primary imbibition characteristics of a reservoir;generating a plurality of constitutive equations based on the pluralityof bounding curves and a plurality of material functions, the pluralityof material functions determined based on material properties associatedwith the reservoir; obtaining a numerical solution of the plurality ofconstitutive equations as capillary-pressure scanning curves, thecapillary-pressure scanning curves representing drainage and imbibitionbehaviors of the reservoir during an oilfield operation; and displayingthe output.
 16. The system of claim 15, wherein the material functionsare pre-determined based on characteristics of the reservoir.
 17. Thesystem of claim 15, wherein the material functions are determined bymatching numerically generated secondary drainage and imbibition curvesto empirical secondary drainage and imbibition curves measured from coresamples from the reservoir.
 18. The system of claim 15, wherein theplurality of constitutive equations comprises:${p_{c} = {{\left( {1 - \alpha} \right)p_{c}^{d}} + {\alpha \; p_{c}^{i}}}},{\alpha \in \left\lbrack {0,1} \right\rbrack},{\frac{\alpha}{S} = {{\eta (\alpha)}\left\lbrack {{{{sgn}\left( \overset{.}{S} \right)}\left( {{\chi\left( {S,\overset{.}{S}} \right)} - \alpha} \right)} + {\gamma\left( {S,\overset{.}{S}} \right)}} \right\rbrack}},{\alpha \in \left\lbrack {0,1} \right\rbrack}$where p_(c) is capillary pressure, p_(c) ^(d) is a primary drainagecurve, p_(c) ^(i) is a primary imbibition curve, α is a hystereticfunction parameter, S is fluid saturation, sgn( ) is sign function, {dotover (S)}=dS/dt is a rate of the fluid saturation, and χ(S,{dot over(S)}), γ(S,{dot over (S)}) and η(α) are at least a portion of theplurality of material functions.
 19. The system of claim 18, wherein theplurality of material functions comprises:${\eta (\alpha)} = \left\{ {{\begin{matrix}{{\alpha\psi}_{0},} & {{{if}\mspace{14mu} \frac{S}{t}} < 0} & {{\text{-}{drainage}},} \\{{\left( {1 - \alpha} \right)\psi_{0}},} & {{{if}\mspace{14mu} \frac{S}{t}} > 0} & {{\text{-}{imbibition}},}\end{matrix}{\chi\left( {S,\overset{.}{S}} \right)}} = \left\{ {\begin{matrix}{{\chi_{1}(S)},} & {{{if}\mspace{14mu} \frac{S}{t}} < 0} \\{{\chi_{2}(S)},} & {{{{if}\mspace{14mu} \frac{S}{t}} > 0},}\end{matrix},{{\gamma\left( {S,\overset{.}{S}} \right)} = \left\{ \begin{matrix}{{\gamma_{1}(S)},} & {{{if}\mspace{14mu} \frac{S}{t}} < 0} \\{{\gamma_{2}(S)},} & {{{if}\mspace{14mu} \frac{S}{t}} > 0}\end{matrix} \right.}} \right.} \right.$ where ψ₀ is a material specificconstant and t represents time.
 20. The system of claim 19, wherein theplurality of material functions further comprises:χ₂(S)+χ₁(S)≡1+(K−1)γ₂(S), γ₁(S)=Kγ ₂(S), where K is a constant.